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Summary: The "Newtonian" theory of spatially unbounded, self-gravitating, pres- 
sureless continua in Lagrangian form is reconsidered. Following a review of the pertinent 
kinematics, we present alternative formulations of the Lagrangian evolution equations 
and establish conditions for the equivalence of the Lagrangian and Eulerian represen- 
tations. We then distinguish open models based on Euclidean space R 3 from closed 
models based (without loss of generality) on a flat torus T 3 . Using a simple averag- 
ing method we show that the spatially averaged variables of an inhomogeneous toroidal 
model form a spatially homogeneous "background" model and that the averages of open 
models, if they exist at all, in general do not obey the dynamical laws of homogeneous 
models. We then specialize to those inhomogeneous toroidal models whose (unique) 
backgrounds have a Hubble flow, and derive Lagrangian evolution equations which gov- 
ern the (conformally rescaled) displacement of the inhomogeneous flow with respect to 
its homogeneous background. Finally, we set up an iteration scheme and prove that the 
resulting equations have unique solutions at any order for given initial data, while for 
open models there exist infinitely many different solutions for given data. 



1 e-mail: ehlers@aei-potsdam.mpg.de 

2 e-mail: buchert@stat.physik.uni-muenchen.de 







1. Introduction 



The Lagrangian theory of gravitational instability of Friedmann-Lemaitre cosmologies 
turned out to be a much more powerful tool for the modeling of inhomogeneities in New- 
tonian cosmology than the standard Eulerian perturbation approach was (for the latter 
see, e.g., Peebles 1980, 1993, and ref. therein). 

Already the general first-order solution of this theory (Buchert 1989, 1992) (which 
contains the widely applied "Zel'dovich approximation" , Zel'dovich 1970, 1973, as a special 
case) has been found to provide an excellent approximation of the density field in the weakly 
non-linear regime (i.e., where the r.m.s. deviations of the Eulerian density contrast field 
6 := p/ph — 1 are of order unity) in contrast to the Eulerian linear theory of gravitational 
instability (Coles et al. 1993, Buchert et al. 1994, Bouchet et al. 1995, Sahni & Coles 1995). 
This appears to be due to the fact that, in contrast to the Eulerian scheme, the Lagrangian 
approximation takes fully into account, at any order, the convective part (v ■ V)v of the 
acceleration and conservation of mass. Another advantage of the Lagrangian equations is 
that they are regular at caustics (where the density blows up), whereas Euler's equations 
break down. Therefore, Lagrangian solutions can be continued accross caustics, i.e., at the 
places where structures form. 

Most recently, the range of application of Lagrangian perturbation solutions for the 
modeling of large-scale structure has been greatly extended by employing filtering tech- 
niques which discard high-frequency modes in the power-spectrum of the initial data, 
and so enable to model highly non-linear stages, even in hierarchical models with much 
small-scale power (Melott et al. 1994, 1995, Weifi et al. 1996). 

In view of these results we think that the power of the Lagrangian description, usually 
applied only to flows under very restrictive conditions (planar, incompressible, etc.), has 
been underestimated. The recent investigation of solutions demonstrates that the com- 
plicated nonlinear partial differential equations which result from the transformation of 
the Eulerian equations to Lagrangian coordinates can be solved in special cases even in 
three dimensions (see Subsection 3.2.3), which was claimed to be impossible in standard 
text books on hydrodynamics discussing the Lagrangian picture. One reason for the pos- 
sibility of constructing solutions lies in the close correspondence of Lagrangian flows and 
classical point mechanics: the Lagrangian coordinates label fluid elements like coordinate 
indices, and in perturbation theory the Lagrangian evolution equations for dust reduce to 
a sequence of ordinary differential equations, as will be shown below. 

For details on the Lagrangian picture of fluid motion in classical hydrodynamics see 
Serrin (1959) and the compilation by Stuart & Tabor (1990). 

We shall treat the initial value problem for the Lagrangian perturbation equations of 
all orders, using a global gauge condition to fix the relation between the background and 
the perturbed flows, and we establish existence and uniqueness of perturbative solutions 
for toroidal (or spatially periodic) models, thus complementing work by Brauer (1992) and 
Brauer et al. (1994). 

Lagrangian perturbation theory has become popular; various authors pursue similar 
studies in relation to the modeling of large-scale structure in the Universe (Moutarde et 
al. 1991, Bouchet et al. 1992, Lachieze-Rey 1993, Gramann 1993, Munshi et al. 1994). 
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For reviews see Bouchet et al. (1995), Bouchet 1996, Sahni & Coles (1995) and Buchert 
(1996a,b). 

Recent efforts concerning the Lagrangian theory in general relativity and in particular 
Langrangian perturbation solutions have been also focussed on evolution equations for fluid 
quantities such as shear and vorticity, the gravitational tidal tensor as the "electric part" 
of the Weyl-tensor, as well as the "magnetic part" of the Weyl-tensor. Supported by the 
classical works by Ehlers (1961), Triimper (1965) and Ellis (1971), a variety of perspectives 
in cosmology have been opened, see the works by Kasai (1992, 1995), Matarrese et al. 
(1993, 1994), Croudace et al. (1994), Salopek et al. (1994), Bertschinger & Jain (1994), 
Bertschinger & Hamilton (1994), Bruni et al. (1995), Kofman & Pogosyan (1995), Lesame 
et al. (1996), Ellis & Dunsby (1996), Bertschinger (1996), Matarrese (1996), Matarrese 
& Terranova (1996), Russ et al. (1996). In these works also the Newtonian limits, or 
analogues, respectively, have been discussed. In a separate note we complement this focus 
by giving a clear-cut definition of the Newtonian limits of the electric and magnetic parts 
of the Weyl-tensor in a 4— dimensional "frame theory" which covers both Newton's and 
Einstein's theory (Ehlers & Buchert 1996). In Newton's theory such fluid quantities are 
expressed in terms of functionals of the trajectories. We emphasize that our point of view of 
a Lagrangian treatment of evolution equations, which was begun with the formulation of a 
closed Lagrangian system for the trajectories by Buchert & Gotz (1987), aims to determine 
fluid quantities explicitly in terms of the trajectory field, and even integrate these quantities 
along the trajectories, if possible, thus, reducing the description to a single dynamical field 
variable. This point of view enables to determine explicitly the evolution of fluid quantities 
without specifying particular solutions for the trajectories. 

The paper is organized as follows: 

In Section 2 we summarize some pertinent facts on the kinematics and dynamics of New- 
tonian self-gravitating flows in the Lagrangian framework. We give an alternative formu- 
lation of the Lagrangian evolution equations in terms of differential forms, we address the 
initial value problem, the problem of existence of solutions, and the equivalence of Eulerian 
and Lagrangian formulations up to the stage when shell-crossing singularities occur. We 
aim to give a self-contained representation of the equations and some additional useful 
relations. Therefore, some equations are reviewed which are not needed in the following 
sections 

In Section 3 we discuss the Lagrangian theory of gravitational instability of the New- 
tonian analogues of Friedmann cosmologies. Here, we give the general perturbation and 
solution schemes at any order and discuss the modeling of space as a 3— torus T 3 as com- 
pared to IR 3 . We give detailed remarks on the interpretation of the perturbation scheme 
and prove uniqueness of the perturbation solutions at any order on the 3— torus. 
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2. The Lagrangian framework 



2.1. Kinematics 

2.1.1. Integral-curves and displacement maps 

Let v[x, t] denote a smooth Eulerian velocity field on IR 3 x [to, t{\. 

We assume that \v\ < V, \dvi/dxj\ < M (indices run from 1 to 3) *. Then there exists a 
unique smooth vector field f(X,t) such that 

±L = v[f{X,t),t] , f(X,t )=:X . (la, b) 

The integral-curves t i— > x(t) = f(X,t) of the velocity field are labelled by the (initial) 

Lagrangian coordinates X; d/dt := d /dt + v ■ V is the total (Lagrangian) time derivative, 
henceforth abbreviated by a dot; a comma (or V) denotes differentiation with respect to 
Eulerian coordinates, and a vertical slash (or Vo) denotes differentiation with respect to 
Lagrangian coordinates; only the latter commutes with the dot. Since dependent variables 
will sometimes be expressed either in terms of Eulerian or in terms of Lagrangian coor- 
dinates, we emphasize the different functional dependence by using the notations [x, t] or 
(X,t), respectively. 

Our assumptions on v imply the following statements (A — G) : 

The integral — curves defined by / do not intersect . (A) 

Since the volume expansion rate 9 := V • v is bounded by 3M, and since (1) gives for the 
Jacobian 

J := det(/ i|fc ) (2a) 



the equation 

we obtain 
Due to (la), 



J(Jt,t) = <£"M*.™ , ( 26 ) 

< e - 3M (*i-*o) < J(X,t) < e 3Wi-to) . ( S ) 

1/1 < v . (C) 



The definition (la,b) of / implies that 

fi } k = v iy ife ]k ; (2c) 



* We employ orthonormal coordinates Xi and use corresponding vector and tensor com- 
ponents; therefore all indices may be written as subindices. 
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therefore, the elements of the deformation gradient Vo/ are bounded, 

|/i|fc| < e 3M (*i-*o) j (jD) 

and 

|/^|<3Me 3M ^- t0 ) . (£) 

— » 

These properties further have the consequences that the displacement map / ( : I h f = 
f(X,t), which sends fluid particles from their initial positions at time to to their positions 
at time t, has the following property: 

ft is an orientation preserving diffeomorphism of IR 3 onto itself , (F) 

(see Appendix A for a proof). 

Let ht denote the inverse of ft', X = h[x,t]. Its Jacobian matrix is given by 

hj,£ = ~^~j€tpq£jrsfp\rfq\s ? (3a) 

and therefore 

\hj,e\ < e 9M ^-*°) . (G) 

So far, we have listed consequences of the definition (la,b) of / in terms of v. Let us 

now, conversely, assume that we have a smooth f(X,t) which has, on IR 3 x [to,ti], the 
properties (A), (£?), (C), (D), (E). Then it is easily established that (F) and (G) also hold, 
and the Eulerian velocity field 

v[x,t] :=f(h[x,t],t) (36) 

is smooth and enjoys boundedness properties of the kind we started with. These remarks 
show under which assumptions the kinematics defined by an Eulerian v[x, t] or a Lagrangian 

f(X,t), respectively, are equivalent; we then call the kinematics regular. 
Remarks: 

(i) The preceding statements remain valid, with some adaptations, if space is modeled not 
as IR 3 , but as a torus T 3 . 

(ii) If, contrary to our assumptions, the velocity field v or / were not bounded, fluid 
particles might escape to infinity in a finite time. If 9 — > — oo sufficiently fast, then J — > 
there, and ft would no longer be locally diffeomorphic; the flow would then develop a 
caustic. If (A) were violated, f t would no longer be injective. In all three cases, (F) would 
fail. 

Under the assumptions discussed above we can also obtain the Eulerian acceleration field 
g = v t t + v • Vv from /: 

g[x,t]:=f(h[x,t],t) . (3c) 
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It is convenient to introduce the following abbreviation: Calculating the Eulerian velocity 
gradient we obtain, with (3a), 

Vi,£ = fi\jhj,£ = -e£ pq J( fi, fp, fq)J~ 1 , (3d) 

where J {A, B, C) abbreviates the functional determinant of any three functions A(X,t), 
B{X,t), C(X,t) with respect to Lagrangian coordinates: 

* A > B >V =:J(A,B,C), 



d(X u X 2 ,X 3 ) 

e.g., for the Jacobian we simply have J = J(f\, f 2 , fz)- 

We now write the curl and the divergence of g in terms of /, using h as a transformation 
from Eulerian to Lagrangian coordinates (hereafter, repeated indices imply summation, 
with k running through the cyclic permutations of 1, 2, 3): 

(V x g) k = epqyJifqJp, f q )J~ 1 , (4a, b, c) 

(V-g) = ^e abc J(f a J b J c )J- 1 . (4d) 

Explicitly, these equations read (summation over j !): 

(V x g)i = J(fj, fi, fj) J -1 , (4a, b, c) 

(V-g)= (j(f u f 2 , f 3 ) + J(f 2 , f 3 , h) + J(f 3 , fi, f 2 j) J' 1 ■ (4d) 

The arguments on the left are x, t, on the right, h[x, t], t. 

Below we give an alternative formulation by using differential forms: 
Let d denote the operator of spatial exterior differentiation acting on functions and forms 
which may be expressed for regular kinematics either in Eulerian (if) or Lagrangian (X) 
coordinates. Then, equations (4) read: 



2 
and 



X (V x g)i e ijk dxj Adx k = g[ij]dxj A dx t = dfiAdfi = d(fidfi) , (4a, 6, c) 



(V • g) dx x A dx 2 A dx 3 = 3df {1 A df 2 A df 3] = d(*fidfi) , {Ad) 



where here * denotes the Hodge star operator with respect to the Euclidean metric dx 2 . 
We shall, however, work with the first form of equation (4d) which turns out to be more 
convenient than the more elegant second form. Also, we shall later use the Hodge star 
operator with respect to the metric dX 2 which coincides with the Euclidean metric dx 2 
only at t = to. The latter operator we shall denote with 
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Recall that the anti-symmetric part taken over 3 indices multiplied by 3 coincides 
with the sum of all cyclic permutations in expressions which involve wedge products, e.g., 

3df {1 A df 2 A df 3] = d'h A df 2 A dh + df 2 A df 3 A dh + df 3 A dh A df 2 . 



2.1.2. Principal invariants of a linear map 

A linear map A : IR 3 — > IR 3 has the following three principal scalar invariants: 

1(A) : = tr(A) , (5a) 
11(A) ■. = \{(tr(Af-tr(A 2 )) , (56) 
1 1 1(A) : = det(A) . (5c) 

For cartesian components, A = (A 1 -) = (Aij). 

In previous work the symbols /, II, and III for the principal invariants of any linear 
map have been used, either with respect to Eulerian or Lagrangian coordinates. The 
kinematical scalars for the expansion, the shear, and the vorticity of the flow f(X, t), which 
we shall use in the present work, can be expressed in terms of the principal invariants (5), 
which we shall do now. 



2.1.3. Relation to kinematical variables 

Let us split the Eulerian velocity gradient (vij) into its symmetric and anti-symmetric 
parts, 

V i,j = V (i,j) + V [i,j] =: ®ij + U ij > ( 6fl ) 

the rate of deformation 9ij and the rate of rotation u>ij. We can split 9ij into its trace-free 
part, the (symmetric) shear tensor a^, and its trace 6, which was introduced already, 

Oij = + -SijO . (66) 

The (anti-symmetric) tensor —u>ij is dual to the angular velocity u, defined as 

uj := x v . (6c) 

— * 

The vorticity tensor uJij = —tijkOOk can be expressed in terms of /, 

UJ ij =V [i,j] = 2 e P q[j^{fi]-.fp-:fq)J~ 1 5 ( 6c 

or, using differential forms, 

uj := —Uijdxi A dxj = dv = d(vjdxj) = dfj A dfj . (6e) 
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The components of uj, u>i = — ^tijkUJjk can be written explicitly as (summation over j !) 

«>i = lj(fjJiJj)J- 1 ■ (6/) 
The magnitudes of shear and rotation are given by 



The preceeding definitions imply 

\ VijVij = u 2 + a 2 + ^9 2 , (7a) 

\ r,.ir hl = -u 2 + a 2 + -9 2 . (76) 

In view of (6) and (7) the principal scalar invariants I, II and III of the tensor (vij) 
are expressible in terms of kinematical scalars, 

I(vij) = Vyi = V • v = 9 , (8a) 

= \ {(vi,if - Vijv jti ) = ^V-(vV-v-v-Vv)=u; 2 -a 2 + ^9 2 , (86) 

11 11 

III(vij) = -VijVj^Vk^ - -K;)Kj^) + g(ui,i) = 3 (viVij)^ 

= - V • ' • v — v • Vv) v + (vV • v — v • Vv) ■ Vv 

= ^-9 3 + 29(a 2 + \u 2 ) + (JijCJjkOki - VijUiUj , (8c) 

where Vij is the matrix with the subdeterminants of Uij as elements. The second equalities 
in (8a-c) show that all invariants can be expressed in terms of divergences of vector fields 
(which has been used and discussed in the context of perturbation solutions - see Buchert 
1994). In obtaining them, the flatness of space is used essentially. 

The velocity gradient Vij = + v [i,j] nas 5 m general, 6 independent scalar invari- 
ants: 

9 , cr , uj , r := -<Jij<JjkO~ki , cfij^i^j , ^ij^jk^i^k , (8d) 

and determines an invariant, ortho normal triad, the eigen-triad of the shear tensor; these 
data together with the 3 Euler-angles of the triad characterize the 9 elements of Vij 
invariantly at any event. 

TruesdelPs invariant, dimensionless vorticity measure (see Serrin 1959) is equal to 

li := : . (8e) 



\A 2 + ¥ 2 



All these kinematical variables can be expressed in terms of / and its derivatives by means 
of eqs. (3). 

It is useful to define the Lagrangian ( "comoving" ) time-derivative of a spatial differ- 
ential form (such as u) in equation (6e)) as the partial t— derivative, taken at fixed Xi, dXi. 
(For the intrinsic, geometrical meaning of this derivative see Appendix B.) 

Then, (6e) implies 

to = dfi A dfi = d(gidxi) = dg . (9) 
Therefore, we have the following kinematical Lemma: 

Let v[x, t] be a (continuously differentiable) velocity field and g = v the corresponding 
acceleration field. Then g is irrotational, V x g = 0, if and only if its vorticity two-form 
u) is conserved in the sense that 

u> = , i.e. , u t = u to . (10) 
(For several equivalent formulations see Appendix B.) 



2.2. Dynamics of self-gravitating "dust" 

So far we considered only kinematical relations which hold for any regular flow field /. We 
now formulate the dynamical equations for Newtonian self-gravitating flows, restricting 
attention to pressureless matter ( "dust" ) throughout this paper. Henceforth the variables 
Xi are to be interpreted as orthonormal coordinates of a dynamically non-rotating frame 
of reference. 

2.2.1. Conservation of mass 

In the Lagrangian framework mass-conservation states that for a regular flow 

g (X,t) = -i r -°Q(X) . (11a) 
J(X, t) 

The Eulerian mass density g can be calculated from (11a) by using the inversion map 
h[x,t]: g[x,t] = g(h[x,t],t). 

Given g(X) > 0, we have shown that under the assumptions of Subsection 2.1.1, g is 
finite and positive for to <t <t\. If, contrary to those assumptions, J — > 0, then g — > oo. 

In terms of differential forms equation (11a) states that the density three-form g d 3 x = 

g dx\ A dxi A dx$ is constant along the flow /: 

gd 3 x = gd 3 X , (116) 
hence -^(g d 3 x) = g d 3 x + g 3dv[i A dx^ A dx 3 ] = (g + gV ■ v)d 3 x = 0, i.e., 

q + q0 = . (11c) 
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2.2.2. Gravitational field equations 

For regular flows, "Newton's" gravitational field equations, generalized by a cosmological 
term, 

Vxj = ; V • g = A - 4nGg , (12a,b,c,d) 

are, in view of equations (4), equivalent to the system of four Lagrangian evolution equa- 
tions (obtained first by Buchert & Gotz 1987 (A = 0) and Buchert 1989 (A ^ 0)): 

J(fjJjJk) =0 , (13a, 6, c) 

(j(/b/2, h) + J(h, fs, fi) + J(fs, h, / 2 )) - A J = -A-kGq . (13d) 
Expressed in terms of differential forms, the Lagrange-Newton system (13) reads: 

dfrAdfr = dtfjdfj) =0 , (13a, b, c) 

and 

3df {1 A df 2 A d/ 3] - A (d/i A df 2 A d/ 3 ) = -4nG°g (dX t A dX 2 A dAT 3 ) . (13d) 

(We keep the numbering (a,b,c) here to remind the reader that these are in fact three 
equations). Equation (13d) can also be written more compactly by using the Hodge star 

operator (with respect to the metric dx 2 ): 

*d(*f j df j ) = A - 4nGg , (13d) 

where g is given by the integral (11a). 

The kinematical Lemma stated at the end of Subsection 2.1.3 shows that in the case 
of "dust", eqs. (12a,b,c) are equivalent to the vorticity conservation law (10) which, in 
this case, acquires the status of a law of gravitational dynamics, dfi A dfi = uj to . In 
particular for irrotational "dust" -flows, u? = 0, the only remaining local law of gravity is 
the divergence law (12d), but the equations dfi A df \ = must not be forgotten! 

The equations (13) are invariant under constant rotations 1Z and time-dependent 
translations T, 

f(X,t)^K-f(X,t) + T(t) , (14a) 

which correspond to the transformations 

x a ' =nt'x b + T a '(t) (146) 

of the Eulerian coordinates. With respect to (14b), the components of the gravitational 
field strength g transform according to 

g a '[x c \t]=n a b 'g b [x c ,t] + f a '(t) . (14c) 

In contrast to the case of isolated systems, where one puts A = and restricts attention 
to inertial frames and Galilean transformations (T = 0), in cosmology the assumption of 
large-scale homogeneity does not allow to single out some coordinate systems as inertial 
ones, and the inhomogeneous term in (14c) unavoidably occurs in transformations relating 
dynamically equivalent coordinate systems (Heckmann & Schiicking 1955, 1956). Then, 
eq. (14c) shows that the gravitational field strength can no longer be considered as a 
spatial vector field independent of the spacetime coordinate system. We shall come back 
to this well-known, but frequently disregarded fact in Subsection 3.1.1. - The arbitrariness 
in the choice of 1Z and T can be restricted or even removed by global conditions depending 
on the solutions considered, as we shall see later. 
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2.2.3. Relations between the Eulerian and the Lagrangian formulations 

The equations (13) are second-order evolution equations for the single dynamical field- 
variable /. An evolution equation for the density is not needed, since g is given explicitly 
by (11). Thus, only three functions of four variables determine the evolution of the system. 
In the Eulerian picture we have seven functions of four variables, e.g., the density, and the 
three components of the velocity and the acceleration field, obeying first-order equations. 

Nevertheless, the regular solutions of the two systems (those with regular kinematics 
in the sense of Subsection 2.1.1) are in one-to-one correspondence, as follows from the 
preceding considerations and has been indicated in (Buchert 1992). More general solutions 
of either system exist, but in general they are no longer equivalent to solutions of the other 
system; see Remark (ii) below. 



Remarks: 

(i) The transition Lagrange — > Euler is simpler than the converse process: in the former 
case, only the equations x = f(X,t) have to be solved "algebraically" for X, whereas in 
the other case, one has to solve the differential equations (1) for /. 

(ii) In writing the first version of equations (13a,b,c,d) we dropped the factor J -1 in front 
of all terms. This is, of course, permitted as long as J ^ 0; it holds in particular for regular 
solutions. Since those equations are regular even at singularities of the system of flow lines, 
i.e., where J = 0, and, in general, J changes sign, one may consider Lagrangian solutions 

o -. 

which have caustics or intersecting trajectories. One may define g[x,t] = J2i ' 

where the sum is performed over all values Xi such that f(Xi,t) = x. Such solutions, 
which contain "multi-dust" regions, are no longer equivalent to Eulerian ones. Their 
physical meaning and validity requires separate considerations and is by no means obvious. 
In particular, they cannot be considered as weak limits of Vlasov-Poisson solutions, since 
in the multi-stream region particles at the same place with different velocities in general 
have different accelerations, which violates the weak principle of equivalence. A general- 
relativistic theory for multi-dust spacetimes which does not suffer from this defect, has 
been outlined by Clarke & O'Donnell (1992). It would seem to be useful to develop a 
corresponding Newtonian theory. Compare also discussions of this problem by Gurevich 
& Zybin (1995). 
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3. Newtonian Cosmology in Lagrangian Form 



3.1. Basic concepts and equations 

3.1.1. Euclidean and toroidal cosmological models 

In Newton's original theory, which was designed and well-defined for isolated systems 
only, as well as in standard versions of "Newtonian" cosmology (see, e.g., Heckmann & 
Schiicking 1955, 1956, or Heckmann 1968), physical space is assumed to be "the" Euclidean 
space based on the manifold IR 3 . For some purposes it is useful or even necessary to model 
3— space as closed, i.e., compact without boundary, as we shall argue in Subsection 3.1.3. 
It is indeed possible to do that without changing any of the local laws so far adopted. 

Since a closed, locally Euclidean 3— space is isometric to the quotient of a flat torus by a 
finite group of isometries* (Kobayashi & Nomizu 1963), we may without loss of generality 
take space to be such a torus T 3 . It is then still possible to cover space at each time 
by finitely many overlapping orthonormal coordinate systems related by transformations 
(14b) with TV 0. 

The inhomogeneous transformation law (14c) for the gravitational field strength can be 
understood by reformulating Newton's theory in covariant spacetime language as initiated 
by Cartan (1923, 1924) and completed by Trautman (1966) (see also the recent work 
on Newton-Cartan cosmology by Rueede & Straumann (1996)). In that reformulation 
the gravitational field is represented as a symmetric, linear connection on spacetime, as 
in General Relativity. It then turns out that there exist non-rotating orthonormal local 
coordinates (t, x a ) such that the only non-vanishing components of the connection are 
given by F^ t . Moreover, the transformations relating these coordinates are those given 
by (14b), and with respect to them the transform exactly like the g a . In fact, the 
free-fall law x a = g a , rewritten as the geodesic equation x a + T% t = 0, shows that we have 
the identity g a = —T^ t , which "explains" the inhomogeneous transformation law and will 
prove useful below. 

3.1.2. Existence of solutions 

Neither the Euler-Newton system, nor the Lagrange-Newton system is a differential system 
to which standard existence theorems apply. The first system is mixed hyperbolic-elliptic, 
while the second is an overdetermined implicit system not fitting into the standard clas- 
sification of PDE theory; the latter may better be considered as an ordinary differential 
equation for the evolution of the time-dependent displacement map. (In this respect, the 
analogous equations of General Relativity are "simpler" (Foures-Bruhat 1958).) Neverthe- 
less, Brauer (1992) succeeded in proving linearization stability of the Euler-Newton system 
at spatially compact (i.e. periodic) Friedmann-like solutions and local-in-time existence 
and uniqueness of solutions which represent finite perturbations of those cosmological mod- 
els, and Brauer et al. (1994) strengthened this result in several ways. The existence and 
uniqueness results established in these papers refer to deviations from a spatially compact 
homogeneous background model which has to be specified, at least partly, for all time 

* in particular, it cannot have the topology of a 3-sphere, a fact which excludes "New- 
tonian" cosmological models based on a 3-sphere 
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and not just by initial data; they do not refer to the total solution (background + per- 
turbation). In fact, "the field equations of the Newton-Cartan theory [a 4-dimensional 
reformulation of "Newton's" theory], unlike the Einstein equations, "are not strong enough 
to determine a solution uniquely in terms of initial data" (Brauer et al. 1994). For this 
and other reasons, work in Newtonian cosmology should be considered as a step towards 
corresponding relativistic considerations. 

Known solutions of the Lagrangian equations include Newtonian analogs of Fried- 
mann's and Bianchi-type general-relativistic cosmological models. Some exact inhomoge- 
neous solutions have also been found (see Subsection 3.2.3). 

3.1.3. Locally isotropic cosmological models 

Those fluid motions which are locally isotropic in the sense that, at any time and for each 
fluid particle V, there exists a neighbourhood on which the field of velocities relative to V 
is invariant under all rotations about V, are characterized by u> = 0, a = 0, V6> = and 
given with our coordinate choice (lb) by 

x = f H (X,t)=a(t)X , a(t ):=l , (15) 

if we conventionally put /j?(0, t) = 0. Such a motion, a Hubble flow, solves the Euler- 
Newton or the Lagrange-Newton system, respectively, if and only if Friedmann's equation 
holds, 

d 2 -e 8ttGq h + A 

= ; e = const. , (lb) 



a 2 



which implies 



a -AnGQ H + h 
a = 3 ' (16) 

where qh = QH(to) a ~ 3 denotes the homogeneous density, and e, A and QH(to) are con- 
stants. Equation (16) holds as well in General Relativity, where the energy constant e is 
related to the Gaussian curvature K at to by e = — K c 2 . Local isotropy implies spatial 
homogeneity, as is well-known. 

Instead of considering the 3-spaces t = const, of the locally isotropic, Friedmann-like 
solutions as globally Euclidean, we may consider the latter as closed, i.e., without loss 
of generality as toroidal, as remarked above. The simplest case arises if we identify all 
those points (particles) whose Lagrangian coordinates differ by integer multiples of some 
constant length L (for the general case see Brauer et al. 1994). In order not to burden 
our equations by powers of L, let us choose L as our unit of length, i.e., put L = 1. 
All particles of such a toroidal universe change their distances in proportion to a(t), the 

locally Euclidean metric is olx 2 = a 2 (t)dX 2 as before, but now the total volume of the 
universe is a 3 (t). Note that this universe is homogeneous and locally, but not globally 
isotropic. The coordinate lines X a = const, correspond to the shortest closed geodesies 
(of length L = 1); geodesies of different directions may be closed and longer, or not closed 
and of infinite length. If we fix an orientation (handedness), the coordinate system (X a ) is 
now intrinsically fixed except for translations and those rotations which map the preferred 
ortho normal triad onto itself. This removes the arbitrariness of 1Z in eq. (14a) except for 
the 9 rotations just mentioned. 
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The toroidal space as a differentiable manifold cannot be covered in a one-to-one, 
bicontinuous manner by a single coordinate system. The coordinates (X a ) used so far 
are coordinates on IR 3 , the covering space of the torus T 3 . In order to see whether the 
gravitational field is well-defined on the spacetime with toroidal space, it is inconvenient 
to use Eulerian coordinates (x a ) and the corresponding g b = — T b t = ^x b ; for then one 
would have to cover T 3 by several overlapping Eulerian coordinate systems and use the 
inhomogeneous transformations to relate the ^-components in the overlap regions. It 
is easier and more elegant to transform the connection components T b t via the geodesic 
equation 'x b — ^x b = to the X "-coordinates. Since x b = a(t)X b , we obtain X b + 2^X b = 
0, for arbitrarily moving test particles (not to be confused with the particles following 
the cosmological flow). Consequently, the non-vanishing components of the gravitational 
connection are T b c = ^5 b . This formula shows immediately that the connection passes 

from IR 3 to T 3 . In fact, instead of working "instrinsically" on T 3 , we may use coordinates 
(X b ) on IR 3 , with the agreement that coordinate values (X a ) differing by integers (iV a ) 
label the same point of T 3 , and provided the relevant fields are periodic. The T b c are not 
only periodic, but translation and rotation invariant due to the homogeneity and local 
isotropy of the model. (This is not obvious in terms of Eulerian components.) 

In Subsections 3.1.5 and 3.2 we shall consider inhomogeneous models as (finite) deviations 
from "Friedmann" -models on T 3 , using "periodic" Lagrangian coordinates {X a ). The 
reason for using T 3 instead of IR 3 is as follows. We shall set up a sequence of perturbation 
equations and show that on T 3 the solutions to these equations to any order exist and 
are uniquely determined by initial data, in accordance with a non-perturbative result of 
Brauer et al. (1994). On IR 3 , however, the corresponding solutions are determined, at 
each order, up to harmonic functions only, i.e., there are infinitely many solutions for the 
same data. 

Uniqueness can also be achieved on IR by restricting the perturbations to be square- 
integrable. Such perturbations, however, contradict large-scale homogeneity. Moreover, 
it is usual to work with periodic perturbations, which can conveniently be represented by 
(discrete) Fourier series. In any case, on T 3 , but not in general on IR 3 , it is possible to 
relate initial and final perturbations unambiguously. 

Remark: 

We can also discuss this problem from a statistical point of view: If one represents the 
typical features of the Universe not by one solution, but by an ensemble, one can maintain 
statistical homogeneity (Bertschinger 1992) in terms of an ensemble consisting of square- 

integrable members, i.e., in terms of perturbations P on IR 3 (introduced below) satisfying 

J d 3 X P 2 (X) < oo. Plancherel's theorem asserts that then the perturbations are also 

A 

square-integrable in Fourier space, i.e., f d 3 k \P\ 2 (k) < oo. Additionally, we may then 
choose the power spectrum of the density perturbations to obey fall-off conditions which 
guarantee square-integrability of the whole random field. Provided that all individual 
members of the statistical ensemble are square-integrable (not merely statistical averages), 

we can set limits on the exponent of a power spectrum of power law form P oc \k\ n : On 
the small-scale end (\k\ — > oo) we have to require n < —3, and on the large-scale end 
(| A; | — > 0), n > —3 (Here we refer to the relations (27a,b) given below and the well-known 
relation between peculiar-velocity and density contrast in the linear regime). Actually, 
the large-scale asymptotics can be satisfied easily, where n ~ +1 according to the COBE 
observations, but the small-scale asymptotics is logarithmically divergent for n = —3, and 
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the maximally allowed slope is n ~ —3 if the spectrum is, e.g., truncated exponentially. 
The latter requirement is at the border of what is allowed in current structure formation 
scenarios. 

Nevertheless, as we have shown in (Buchert & Ehlers 1996), spatially closed universes 
(i.e., those which are compact without boundary) are singled out as the only generic models 
in which the averaged variables of inhomogeneous fields represent homogeneous solutions. 
Thus, the toroidal universe is the simplest among those Newtonian cosmologies. 

3.1.4. Average properties of general 

inhomogeneous cosmological models 

Following Buchert & Ehlers (1996) we discuss spatial averages of inhomogeneous New- 
tonian cosmological models by deriving the general expansion law which is obtained by 
averaging Raychaudhuri's equation (Raychaudhuri 1955): 

9 = A - 4ttGq - \e 2 + 2(uj 2 - a 2 ) . (17) 
(Differentiation of the expansion scalar 9 with respect to the time yields 

= Vi,i,t + VjVi,i,j = v ijt ,i + (vi,jVj),i ~ Vi,jVj,i = 9i,i + 2 ^ 2 - 2a 2 - ^6> 2 . (17') 

In view of (12d) we obtain (17).) 

Equation (34) shows that if, on one trajectory, ^A + u 2 < 2ttGq + a 2 (in particular, 
if A = and u = 0) and 9(t') ^ 0, then there exists an instant of time t" such that 
sgn(*' -t") = sgn(0(£')), \f - t"\ < ^;linw^(t) = ttm t _+ t ,,\0{t)\ = oo. 

Let us consider an arbitrary "comoving" (Lagrangian) volume V(t) =: a|,(t) of a spatially 
compact portion V(t) of the fluid; it changes according to 

v= ^- [ ofx = [ d 3 x j = [ d 3 x e , 

which may be written 

«>)v = I = 3^ . (18) 
V ax> 

Here and in the sequel, (A)x> = y Jx> d 3 % A denotes the spatial average of a (spatial) tensor 
field A on the domain V{t) occupied by the amount of fluid considered, and op is the scale 
factor of that domain. 

The average of Raychaudhuri's equation may then be written (Buchert & Ehlers 1996): 
3^ + 4ttG^ - A = \ {{9 2 ) v - (9)1) +2(u 2 - a 2 ) v . (19) 

(ID CLjy o 

We have used the definitions (6g,h). Equation (19) shows that the presence of inho- 
mogeneities affects the expansion law which only coincides with Friedmann's law (16'), 
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ax> = a, provided shear, vorticity and fluctuations of the expansion scalar vanish or cancel 
each other, respectively. 
Introducing the averages 

:= (9)t> ; Sij := (o-ij)v ; &ij ■= (uij)v , (20a, 6, c) 

we define a linear "background velocity field" V on V by Vi := HijXj with 

Vi,j = + ^QSij + Qij =: Hij . (20d) 

(Note that all average variables, like a(t), 0(£), and dij(t), depend on content, shape 

and position of the spatial domain T>.) 

While the velocity fields v and V depend on the choice of a non-rotating frame of 
reference (c/. eq. (14b)) and are consequently not global vector fields on a toroidal model, 

the peculiar velocity field, defined as u := v — V, always is a global vector field. Splitting 
expansion, shear and vorticity into their (time-dependent) average parts and deviations 
thereof, 

9 = + 9; Uij = + &ij ; Uij = Q,^ + (% , (21a, 6, c) 

equation (19) can be cast into the form 

3 ^ + AnG M_ _ A = 2(0 2 - E 2 ) + -{9 2 ) v + 2{u 2 - & 2 ) v . (22) 
av of, 3 

(The averages (9)x>, (^ij)v and {uJij)v vanish by definition.) 
Using (8b) for the peculiar-velocity gradient (uij), 

2 - 9 2 + 2{u 2 - a 2 ) = V • [u(V -u)-(u- W)u] , 
we finally arrive at the remarkably simple general expansion law: 

3^ + 4nG^- - A = 2(0 2 - E 2 ) + (V • [u(V ■ u) - (u ■ V)u])v ■ (23) 

The last term in (23) is, via Gaufi's theorem, a surface integral over the boundary of V. 
In case of a toroidal model we may choose T> to be the whole torus. Thus, on the torus, 
we obtain the global expansion law (in agreement with the result of Brauer et al. (1994)): 

3— +4ttG^-A = 2(0 2 -E 2 ) ; V = T 3 . (23') 

(I'D O'j) 

This law, combined with the linearity of the velocity field V, can be used to determine all 
homogeneous, in general anisotropic Newtonian models either on IR 3 or on T 3 , in Eulerian 
or Lagrangian form (for models on IR 3 in Eulerian form, see Heckmann & Schiicking 
(1959)). 
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The point of this subsection was to show how these models arise by spatially averaging 
arbitrary inhomogeneous models, provided either space is compact or, if for V —>■ IR 3 , the 
last term in (23) vanishes. 

In the remainder of this paper we restrict ourselves to models having locally isotropic 
backgrounds, i.e., where Ejj = fijj = 0; then, the average motion is a Hubble flow whose 
expansion is described by Friedmann's law (16'). 

3.1.5. Inhomogeneous cosmological models 

as deviations from locally isotropic ones 

We wish to consider periodic or toroidal inhomogeneous models which are isotropic (and 
hence irrotational) on average on some large scale. As shown in the last subsection, the 
requirement of periodicity implies that the spatially averaged density 

< e > > TS (() : = kf^l = M^L = ^L. (24 ) 
f^cPX J(X,t) V(t) o?(t) 

of any such model is related to a(t) by Friedmann's equation (16) with some constants 
e, A, QH(to) (which are then uniquely determined). Thus, we can associate with any 
inhomogeneous model its toroidal locally isotropic background model defined by qh '■= {q)t 3 
and a(t) via eqs. (15), (16), as described in Subsection 3.1.3. 

To describe inhomogeneous cosmological models we define the deviation p of the dis- 
placement map / of the inhomogeneous model from the background model fu- 

f = f H +p(X,t) ; p(X,t ):=6 . (25a, b) 

It is convenient to introduce periodic rescaled Eulerian coordinates *, q := x/a(t) and the 
corresponding deformation field F, q = F(X,t); F(X,t ) = X. Then, the equations (25) 
read: 

F = X + P(X,t) ; P(X,t ):=6 , (26a, b) 

where P = p/a(t). P t : IR 3 — > IR 3 is periodic and may be interpreted as the (conformally 
rescaled) displacement of the particles of the perturbed flow relative to those of the unper- 
turbed flow. It is considered the fundamental object of Lagrangian perturbation theory 
hereafter. 

To fix the (fictitious) mean displacement of the perturbed flow relative to the un- 
perturbed one ("identification gauge condition"), we require, without loss of generality, 
besides (26b) for all t: 

[ d 3 X P(X, t) = . (26c) 

It fixes the choice of T in equation (14a) and is essential for the uniqueness of Newtonian 
solutions, as we shall see later. Note that (26c) can also be written (g/g P)j3 = so that, 
if g is nearly constant, {g P)j3 0, a center-of-mass condition. 



i.e., Lagrangian coordinates of the background flow 
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The displacement vector P determines the peculiar-velocity u and the peculiar- 
acceleration w by: 

u := v x = aP ; u = P(to) , (27a) 

a 

w ■= g- -x = aP + 2aP ; k = P(t ) + 2a(t )P(t ) , (276) 
ct 

o o 

where u and w are the initial data for peculiar-velocity and peculiar-acceleration, respec- 
tively. (Note that while P, u, w are global vector fields on T 3 , the Hubble velocity ^x and 
v are defined only locally with respect to some "origin".) 

Below we shall use the corresponding one-forms denoted by U = UidXi and W = WidXi, 
and for the time-dependent perturbation P = PjdXj. 

Let us now write down the equations which the displacement P has to obey. Inserting 
(26a) into the once integrated Lagrangian evolution equations (13a,b,c) results in 

dPi A (dX l + dPi) = a~ 2 & = d(a~ 2 U) . (28a, 6, c) 

The latter equality follows from (6e) and the fact that the Hubble-velocity is assumed to 
be irrotational. The last equation may be rewritten as 

d{p + PidPi - a" 2 u| = . (28a, 6, c) 

Note that there is no cubic term in these equations. 

Inserting (26a) into (13d), and defining the operator V := jj^z + 2H-^ and the function 
b := 3^ - A, we obtain 

b dX 1 A dX 2 A dX 3 + (V + b)SdP {1 A dX 2 A dX 3] + (V + 2b)3dP [1 A dP 2 A dX 3] 

+ (-V + 6)3dP[i A dP 2 A dP 3] = — ^g rfXi A dX 2 A dX 3 . (28a") 

(Remember that expressions of the form 3dA[i A dA 2 A d^4 3 ] are equal to the sum of all 
cyclic permutations: J2ijk dA-i A dAj A dAk.) 

— * 

Since this equation holds for the background, P = 0, the terms independent of P 
cancel, and we are left with the equation 

(V + 6)3dP[i A dX 2 A dX 3] + (V + 2b)3dP {1 A dP 2 A dX 3] 
+ {\v + b)3dP {1 A dP 2 A dP 3] = ~ 4lT f 5 ° Q dX 1 A dX 2 A dX 3 , (28d) 

where Sg = g — ^ is the (finite) initial deviation from the homogeneous density gn = 
Q°Ha~ 3 ; f T3 d 3 X5°g = 0. 

In what follows we shall use the Hodge star operator with respect to the metric dX 2 . 
Therefore, we indicate it with a big star (*) to avoid confusion with the Hodge star 
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operator used in previous equations. (The following identities are useful: %-d 3 X = 1, 
(*)2 = l, d*d* = *d*d = A .) 

Operating with * on (28d) and using AttGSq = %-d^W ', gives 



*d|(£> + &)*P + (V+2b)3P [1 AdP 2 AdX 3] + (-V+b)3P [1 AdP 2 AdP 3] -a- 3 *wj =0 . 

(28d) 

Here, the linear term is purely longitudinal. 

The equations (28a,b,c,d) with the initial conditions (26b) govern inhomogeneous models. 
In more familiar vector notation the equations (28a,b,c,d) have the form: 

^(V xP) = F(dP t , dPj) + a" 2 V x °u ; 

(V + 6)(V • P) = Q{Pi, dPj, Pi, Pi) + a" 3 V • . 

The r.h.s.'s contain no terms linear in P or its derivatives, and they contain no derivatives 
with respect to t or Xi of higher order than on the l.h.s. . Therefore, these equations 
lend themselves to solution by iteration. For that purpose, the condensed differential form 
notation is more convenient than vector notation, however. 



3.2. Lagrangian perturbation theory 
3.2.1. The perturbation scheme 

Since we have only one dynamical object in the problem (the one-form P), a Lagrangian 
perturbation scheme on Friedmann-Lemaitre backgrounds can be set up by inserting into 
eqs. (28) for P a formal power series, 

oo 

P = £ m P (m) , (29) 

m=l 

to obtain a sequence of equations for the p( m ) at order m. We thus obtain the following 
system of 4m equations: 
For m = 1 we have 

dpW = d(a" 2 UW) ; (30a, 6, c; to = 1) 



For to > 1 we have 



d*{[p + 6]p (1) } = d(a- 3 *wW) . (30d;TO=l) 

rf |p(m)| =rfT (m) . (30a, 6, c; to > 1) 

d*{[p + &]p (m) } = d*S (m) . (30d;m>l) 
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The 2m source terms (one-forms) and T*™) can be read off eqs. (28). They depend 
on P ( ^; I < m: 

m — l 

T (m) = _J2 pWdP^-^i + a" 2 U (m) , (31a; m > 1) 

t=i 

m — l 

^ S (m) = _ (V + 2o)3P (£) [1 dP (rW) 2 A dX 3] - 
1=1 

V (^P + 6)3P W [idP (p) 2 AdP (9) 3] +a" 3 *W (m ) . (316;m>l) 

Starting at the third order, the source terms contain products of perturbation solutions of 
different orders, (compare Buchert 1994, eqs. (4)). 

3.2.2. General solution scheme 

To solve the equations (30) with the source terms (31), we decompose the p( m )'s as well 
as the initial values U and W non-locally into their longitudinal and transverse parts (see 
Appendix C), 

p(m) = p(m) L + p (m) T ? ^2a) 
Jj(m) = jj(m) L + jj(m) T ^ (33^ 
W (m) = W (m) L ^ ( 32c ^ 

taking into account that the harmonic parts vanish because of the gauge condition (26c) 
and eqs. (27), and remembering that dW = 0. 

We prescribe, without loss of generality, that the initial density perturbation and thus 
W be of first order, 

6q = 6qM =:q h 6 ; = W , (33a, b) 

o 

where 5g denotes the initial density perturbation, and 5 the initial (conventional) density 
contrast. 

Equation (26b) requires, for all m, 

P (m) (X,t ) :=0 . (33c) 
Finally we require, also without loss of generality, 

P(X, t ) = P (1) (X, t ) = U(X) . (33d) 

The unique solutions of the perturbation equations having these initial data are obtained 
as follows. 
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Equations (30a, 6, c; m = 1) say that 



A:=pW T -a- 2 UW T 
is both closed, dA = 0, and co-exact, hence it vanishes (see Appendix C); therefore 

P (1)T (Xi) = u T (i) f4L • ( 34 ^ M 

Jto a K 1 ) 

Eq. (30d; m = 1) similarly implies 

(V + b)P (1)L (X,t) =a- 3 W(X) . (34d) 

The solution to this ordinary differential equation obeying the initial conditions (33) is 
uniquely determined by the data W(X) and XJ L (X). 

For m > 1 we obtain from (30a,b,c): 

P^ T (X,t)= [ dt'T^ T (X,t') ; (35a, 6, c) 

and from (30d): 

(V + b)P {m)L (X,t) = S (m)L . (35d) 

The solutions to eqs. (35) are uniquely determined by their sources (31), since they are 
required to have vanishing initial values. 

Remarks: 

(i) The solutions at any order m are well-defined and unique on E x T 3 as long as the 
background is free of singularities. In general they will develop "multi-dust" regions. 

(ii) The solutions at any order m separate with respect to Lagrangian coordinates X and 

time t; p( m )(X,t) = ^ a A^ n \x)B^ n \t). This property follows from the structure of 
the perturbation scheme, since the first-order solutions separate and, at each step, only 
linear ordinary differential equations with respect to t have to be solved. The time- 
dependent coefficients are determined solely by the background, while the X-dependent 
factors depend on the initial data. 

(iii) The first-order solution depends locally on the data U and W in the sense that the 
factors A£\x) at some value X depend only on U and W at the same X. On the other 

o 

hand, W depends non-locally, via a solution of Poisson's equation, on 5. Each further step 
involves the determination of T( m ) and S^ m ^ from T^ m ^ and S^ m \ respectively, which 
again requires to solve Poisson equations. Thus, the X-dependent factors in p( m ) depend 
non-locally on the data U and W for m > 1. The trajectory of each "dust particle" at 
any order of approximation depends globally on the initial data, even at times close to the 
initial time, just as in Newtonian dynamics of systems of finitely many particles. This is 
in contrast to General Relativity, where the evolved fields at some spacetime point depend 
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only on the initial data within the causal past of that point. (For GR "dust" solutions this 
has first been shown by Foures-Bruhat 1958.) 

(iv) Since all relevant functions are defined on T 3 , they can be represented by discrete 
Fourier series. Since the sources for the higher-order terms are products of lower-order 
ones, the higher-order terms will change on smaller spatial scales than the lower-order ones, 
and their time-dependent factors will contain (positive and negative) powers of those of 
the first-order solution which generates the higher-order ones. 

(v) If the perturbation scheme is applied to fields on IR 3 rather than on T 3 , at each step a 
harmonic contribution to p( m ) has to be chosen arbitrarily. (This is due, of course, to the 
form of eqs. (12)). Then, there are infinitely many perturbative solutions for given initial 
data; hence, it makes no sense to ask which fields evolve from which data. 

(vi) The equations (34) suggest that it is convenient to introduce a new time-variable T 
(taken to be dimensionless): 

dT := — — . (36a) 

to a 2 (t) 1 ' 

This variable has been very useful for the purpose of finding solutions for "non-flat" 
backgrounds (see: Shandarin (1980), Buchert (1989, Appendix A), Bouchet et al. 1995, 
Catelan 1995). With this time-variable solutions of (16) for A = have the simple form: 

a(T) = « . (366) 
Also the time-dependent operator in front of the longitudinal part simplifies (A ^ here): 

t 2 (V + b) = -^-4nGp H a . (36c) 

(Compare: Buchert (1989, Appendix A) for the Lagrangian equations as well as all relevant 
cosmological variables and parameters expressed in terms of T) . 



3.2.3. Explicit solutions 

(Not in chronological order of their derivation.) 

Known solutions comprise the general first-order solution (Buchert 1992) for an "Einstein- 
de Sitter" background, which includes rotational flows and the "Zel'dovich Approximation" 
(Zel'dovich 1970, 1973) as the special case U T = 0, U L = Wt . 

For irrotational flows the solution for all backgrounds with A = can be found in 
(Buchert 1989) including generalizations of Zel'dovich's approximation obtained by Shan- 
darin (1980). 

For most of the background solutions including a cosmological constant, closed-form 
expressions are given in (Bildhauer et al. 1992), where a general procedure to obtain the 
"Zel'dovich Approximation" for all backgrounds is outlined. 

Interestingly, for restricted initial data, the first-order solutions turn out to be exact 
three-dimensional solutions (Buchert 1989) including the general plane-symmetric solution 
given earlier by Zentsova & Chernin (1980). These solutions contain caustics. (For related 
exact solutions see Buchert & Gotz (1987), Barrow & Gotz (1989) and Silbergleit (1995).) 
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At second order all irrotational solutions on an Einstein-de Sitter background are 
known for initial data which admit a functional dependence of initial peculiar-velocity 
and peculiar-gravitational potentials (Buchert & Ehlers 1993). A subclass of these solu- 
tions for the special case U T = 0, U L = Wto is discussed in Buchert (1993). For the 
same initial data the third-order solution on an Einstein-de Sitter background is given by 
Buchert (1994), the fourth-order solution by Vanselow (1995); see Sahni & Coles (1995) 
and Buchert (1996a,b) for reviews. 

Lagrangian perturbation solutions and their applications have also been derived and 
applied by Bouchet & collaborators (for a review see Bouchet et al. (1995), where references 
to solutions with "non-parabolic" cosmological backgrounds at second (Bouchet et al. 
1992) and third order for the leading time coefficient (the particular solutions) can be 
found). Moutarde et al. (1991) gave a third-order approximation on an Einstein-de 
Sitter background for special symmetric initial data. For these data a (slightly different) 
solution has been derived from the generic solution by Buchert et al. (1996). The general 
irrotational second-order solution for "non-parabolic" cosmological backgrounds with zero 
cosmological constant has been derived by Vanselow (1995). Also Munshi et al. (1994) 
discuss the leading terms of the third-order solution of Buchert (1994), and Catelan (1995) 
derives and discusses the third-order solution for "non-parabolic" backgrounds. 

The main difference between most of these works and our approach is that we con- 
sistently work within the Lagrangian framework, i.e., we express all equations in terms of 
the single dynamical field / before solving them. Hence, we avoid mixing of Lagrangian 

— * 

and Eulerian representations. The only perturbed field is / in Lagrangian space; all Eu- 
lerian fields are calculated therefrom. The velocity field is determined perturbatively, the 
corresponding mass and the vorticity is exactly conserved in our perturbation solutions. 

The fundamental question whether these perturbation solutions converge to or, at 
least, approximate exact solutions remains open. 
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APPENDIX A 



Under the assumptions stated at the beginnning of Subsection 2.1.1, the map f t : D 3 — > D 3 ; 
D G {K, T}; t fixed; (to < t < ti) is a diffeomorphism. We first show that / t is injective, 
and then that it is surjective. Since f t is a local diffeomorphism because of J > 0, this 
establishes the claim. 

Injectivity follows immediately from the fact that different integral-curves of a vector field 
are disjoint. 

To establish surjectivity we notice the following: 

Since f t is a local diffeomorphism, the image /t(D 3 ) is open. It is also closed; for let 

Xi = ft(Xi) be a sequence of images which converges to xq, Xi — > xq. Then, the set {Xi} 
is bounded since {xi} is, and distances can change during [to,t] at most by 2V\t — to\. 
Therefore, a subsequence of {Xi} converges to some point Xq. Continuity of ft then 
implies that x = f t (X ) e /t(© 3 )- Thus, / t (D 3 ) is both open and closed in D 3 , hence 
equal to D 3 . 



APPENDIX B 



We here give an invariant meaning to the "time-differentiation" of differential forms which 
was used in the main text (the reader may consult standard textbooks on differential forms, 
e.g., Schutz 1980), and we collect different versions of the vorticity conservation law u> = 0. 

Lie-derivative 

We defined the operator ' on spatial differential forms as partial differentiation with 
respect to t for fixed X. In Newtonian spacetime 1R x IR 3 or IR x T 3 , a velocity field v[x, t] 
determines a world velocity field, 

d d d , . 

dt = di + V *dx-- (BA) 



If we use Lagrangian coordinates (X, t) on spacetime, the vector field ^ has components 

(0,1). Therefore, in these coordinates, Lie-differentiation with respect to ^ amounts to 
partial differentiation with respect to t. This shows that 

L_d_A = A (B.2) 



dt 



for all "spatial" differential forms, i.e., differential forms not containing dt, and gives the in- 
variant meaning of ' . This time-derivative commutes with spatial exterior differentiation, 
d. 

We now list some equivalent versions of the vorticity conservation law 

us = , (B.3a) 
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since different versions appear in the literature and are useful for different purposes (for 
all these relations it is necessary that the force is conservative, i.e. the gravitational field 
strength g is irrotational). 

The vector form of (B.3a) reads: 

u = u • Vv — uJV • v . (B.3b) 

We can integrate uj along the integral-curves / to obtain Cauchy's integral (see, e.g., Serrin 
1959, Buchert 1992), 

u = (O-V f)J- 1 ■ (5.3c) 

Equation (B.3c) shows that the vorticity blows up at points of (formally) infinite density 
(J = 0) for generic initial data (see Buchert 1992 for a proof). This implies that caustics 
are associated with strong vortex flows in their vicinity (see also the detailed discussion by 
Barrow & Saich 1993). 

In terms of kinematical variables, the vorticity law reads: 

2 

uji = —-Oui + TijUj . (B3.d) 



APPENDIX C 



In order to make this paper self-contained and to fix our notation we here collect some 
well-known facts about decompositions of vector fields on IR 3 and T 3 , respectively, both 

furnished with the standard flat (Lagrangian) metric dX 2 . 

On IR , any smooth vector field P can be decomposed into a gradient (longitudinal) 
part and a curl (transverse) part, 

P = P L + P T = V t/ + V x A , V -A = . (C.l) 

— * 

Such a decomposition always exists, whether or not P falls off at infinity; but it is not 
unique: if if is a harmonic field, i.e., a field satisfying Vo • H = and Vo x H = 0, then 

P = (V U + H) + (V x A- H) 

gives another representation of the type (C.l), since H = VoV> = Vo x £?, and in this way 
all such representations are obtained. If P as well as the parts P and P are required to 
be square integrable (e £ 2 ), i.e., f d 3 XP 2 < oo, the decomposition (C.l) is unique; square 
integrable harmonic fields do not exist on IR 3 (Dodziuk 1979). Then one can speak of the 

longitudinal, or the transverse part of P, respectively. 
On T 3 , one has a unique decomposition: 

p = V U + V x A + H , (C.2) 
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where the harmonic part H is constant on T 3 (see the remark below) and given by: 

H= I d 3 X P . (C.3) 

— * 

The potentials U and A can also be fixed uniquely by requiring: 

I d 3 X U = , I d 3 X A = , V • A = . (C.4) 
JT 3 </t 3 

Note that, on T , being longitudinal means not only that Vo x P = 0, but in addition that 
/ d 3 X P = 0. Similarly, transversality requires Vo • P = and vanishing average. 

It is convenient to re-express these facts in the language of differential forms rather 
than that of vector fields. Writing P = PidXi for the one-form (covector) associated with 

— * 

P, the form-analogs are: 

P = P L + P T + P H = dU+ *dA + H , (C.2') 

where A and H are one-forms, the longitudinal part is an exact form, the transverse part 
a co-exact form, and the harmonic part a harmonic form, which is determined by 

I d 3 X P = H , (C.3') 

and one may impose 

/ d 3 X U = , / d 3 X A = , d*A = , (C.4') 

JJ3 JT 3 

where in all equations * denotes the Hodge star operator with respect to the metric dX 2 . 

The integration of the perturbation equations in Subsection 3.1.4 is based on the 
following two facts: If a co-exact form P T is closed, dP T = 0, it is the zero-form, P T = 0. 
If an exact form P L is co-closed, d^P L = 0, it is the zero-form, P L = 0. These facts 

follow from the foregoing statements and equations. 
We also recall that Poisson's equation, 

A U = 4nGg , (C.5) 

is soluble on T 3 if and only if J J3 d 3 X g = 0. The solution is then unique except for an 
additive constant which may be fixed by demanding: 

/ d 3 X U = . (C.6) 

Jt 3 

For proofs see, e.g., Warner (1971). 
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Remark: The only harmonic vector-fields H are the constant ones. To see this, we recall 
the vector-identity 

A H = V x (V x H) - Vo(V • H) . (C.7) 

It shows that a harmonic vector field obeys Laplace's equation. Then, its components 
Hi (i = 1, 2, 3) are harmonic functions. For each component, we can apply Green's formula, 

/ H i A H i = [ HiVo(VoHi)= [ {V (^V ^) - (V ^) 2 } 

JJ3 JJ3 Jj3 

g T 3 OH Jj3 

Since the scalars Hi are harmonic, the left-hand-side of the identity (C.8) vanishes. Since 
the torus T 3 has no boundary, we finally conclude 

/ (V #0 2 = , (C.9) 
</t 3 

or, VgHi = 0. Hence, Hi = const.. 
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